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Abstract 

In this paper we present a numerical scheme for the resolution of matrix Ric- 
cati equation, usualy used in control problems. The scheme is unconditionnaly 
stable and the solution is definite positive at each time step of the resolution. 
We prove the convergence in the scalar case and present several numerical ex- 
periments for classical test cases. 
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1) Introduction. 

• We study the optimal control of a differential linear system 

(1.1) ^=Ay + Bv, 

where the state variable ?/(•) belongs to IR"^ and the control variable v{») 
belongs to IR"^ , with n and m be given integers : 

(1.2) ?/(t)eIR", v{t)eJR'^. 

Matrix A is composed by n lines and n columns and matrix B contains n lines 
and m columns. Both matrices A and B are fixed relatively to the evolution in 
time. Ordinary differential equation (1.1) is associated with an initial condition 

(1.3) y{0) = yo 

with yo be given in H"^. Morever the solution of system (1.1) (1-3) is parame- 
trized by the function v{») and instead of the short notation y{t), we can set 
more precisely 

(1.4) y{t) = y{t; yo,v{*)) . 

The control problem consists of finding the minimum u{t) of some quadratic 
functional J(») : 

(1.5) J{u{.)) < J{v{.)), yv{.). 

The functional J(») depends on the control variable function v{»), is additive 
relatively to the time and represents the coast function. We set classically : 

(1.6) J{v{.)) = {Qy{t)Mmt+\j^ {Rv{t)Mmt+\{Dy{ny{T)). 

Functional J(«) is parametrized by the horizon T > , the symmetric semi- 
definite positives n by n constant matrices Q and D : 

(1-7) iQy,y) > 0, V2/€IR^ y^O, 

(1.8) {Dy,y)>0, V?/ € IR", 2/^0. 

and the symmetric definite positive m by m constant matrix R : 

(1.9) {Ru,u)>0, \fueJR'^, uj^O. 

• Problem (1.1) (1-3) (1.5) (1.6) is a classical mathematical modeling of linear 
quadratic loops for dynamical systems in automatics (see e.g. Athans, Falb 
[AF66], Athans, Levine and Levis [ALL67], Kawakernaak-Sivan [KS72], Faurre 
Robin [FR84], Lewis [Le86]). When control function v{») is supposed to be 
squarely integrable {v{») e L^(]0, T[, IR"^)) then the control problem (1.1) (1.3) 
(1.5) (1.6) has a unique solution u{t) (see for instance Lions [Li68]) : 
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(1.10) u{t) €L2(]0,T[, K^). 

When there is no constraint on the control variable the minimum of functional 
J(») is characterized by the condition : 

(1.11) dJ{u).w = 0, V«;€L2(]0,T[, IR"^), 

which is not obvious to derive because y{») is a function of v{»). We introduce 
the differential equation (1.1) as a constraint between ?/(•) and v{») with the 
associated Lagrange multiplayer p . We set : 

(1.12) C{y,v-p) = J{v) - ^ L'^^ - Ay- Bvj dt 

and the variation of £(•) under an infinitesimal variation Sy{») , 6v{») and 
Sp{») of the other variables can be conducted as follow : 

6C= [ iQy,dy)dt+ [ {Rv,6v)dt + {Dy{T),dy{T)) 
Jo ^ Jo ^ 

^ {6p,^-Ay-Bv)dt - ^ (P^^ - - ^Hdt 



[ {Qy,Sy)dt + [ iRv,6v)dt + (Dy{T),6y{T)) 
Jo Jo 



(1.13) 



/ (A^p,6y)dt + / {B^p,6v)dt, and 
Jo Jo 

5C = [ + Qy + A^p, 5y) dt + f {Rv + B^p,5v) dt 
Jo dt Jq 

- J (Sp, ^ - Ay - Bv) dt + (Dy{T) - p{T),6y{T)) 



because Sy{0) = when the initial condition (1.3) is always satisfied by the 
function (y + 6y){») . 

• The research of a minimum for J(») (condition (1.11)) can be rewritten 
under the form of research of a saddle point for Lagrangien C and we deduce 
from (1.13) the evolution equation for the adjoint variable : 

(1.14) ^ + A^p + Qy = 0, 

the final condition when t = T ^ 

(1.15) p{T) = Dy{T) 
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and the optimal command in terms of the adjoint state />(•) : 

(1.16) Ru{t) + B^p{t) = 0. 

We observe that the differential system (1.1) (1.14) joined with the initial condi- 
tion (1.3) and the final condition (1.15) is coupled through the initial optimality 
condition (1.16). In practice, we need a linear feedback function of the state 
variable y{t) instead of the adjoint variable p{t) . Because adjoint state p{») 
linearily depends on state variable y{») we set classically : 

(1.17) p{T) = X{T-t)y{t), 0<t<T, 

with a symmetric n by n matrix X{») , positive definite for t > (see e.g. 
[KS72] or [Le86]). 

(1.18) ^it) is a symmetric nxn definite positive matrix, t>0. 

The final condition (1.15) is realised for each value y{T) , then we have the 
following condition : 

(1.19) X{0) = D, 

and introducing the representation (1.17) in the differential equation (1.14) and 
(1.1) we obtain : 

-^{T-t)y{t) + X{T -t)[Ay{t) + Bu{t)] + X{T -t) y{t) ^ Qy{t) = 0. 

We replace the control u{t) by its value obtained in relation (1.16) and we 
deduce : 

'^^■{T-t)+X{T-t)[A + B{-R-^)B^X{T-t)]+A^X{T-t)+Qy{t) = 0. 



dt 

This last equation is realised for each state value y{t) . Replacing t by T — t 
in this equation, we get : 

(1.20) (XA + A^X) + XBR-^BX - Q = 0, 

which defines the Riccati equation associated with the control problem (1.1) 
(1.3) (1.5) (1.6). 



• In this paper we study the numerical approximation of differential system 
(1.19) (1.20). Recall that datum matrices Q, D and K, with K defined ac- 
cording to : 

(1.21) K = BR-^B^, 

are nxn symmetric matrices, with Q and D semi-definite positive and K 
positive definite ; datum matrix ^ is an n by n matrix without any other 
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condition and the unknown matrix X{t) is symmetric. We have the following 
property (see e.g. [Le86]). 



Proposition 1. Positive definitness of the solution of Riccati equation. 

Let K, Q, D,A be given n x n matrices with K, Q, D symmetric matrices, 
Q and D positive matrices and K a definite positive matrix. Let X{») be the 
solution of the Riccati differential equation : 



(L22) (XA + A^X) + XKX - Q = 

kXL 

with initial condition (L19). Then X{t) is well defined for each t > , is 
symmetric and for each t > , X{t) is definite positive and tends to a definite 
positive matrix X^o as t tends to infinity : 

(1.23) X{t) — > X^ if t — ^ oo . 

Matrix X^o is the unique positive symmetric matrix which is solution of the 
so-called algebraic Riccati equation : 

(1.24) -{XA + A^X) + XKX - Q = 0. 

• As a consequence of this proposition it is usefull to simplify the feedback 
command law (1.16) by the associated limit command obtained by taking t — > 
oo , that is : 

(1.25) v{t) = -R-^B^Xooy{t), 

and the differential system (1.1) (1-25) is asymptotically stable (see e.g. [Le86]). 
The pratical computation of matrix X^o with direct methods is not obvious 
and we refer e.g to [La79] for a description of the state of the art. If we wish to 
compute directly a numerical solution of instationnary Riccati equation (1.22), 
classical methods for ordinary differential equations like e.g the forward Euler 
method : 

(1.26) ^(^.+1 - Xj) + XjKXj - {A^X, + X,A) - Q = 0, 

or Runge Kutta method as we will see in what follows fail to maintain positivity 
of the iterate Xj_^i at the order (j ' + 1) : 

(1.27) {Xj+ix,x)>0, yxeW, Xy^O, 

even if Xj is positive definite if time step At > is not small enough, see e.g. 
[Sa97]. Morever, the stability constraint (1.27) is not classical and there is at 
our knowledge no simple way to determine a priori if time step At is compatible 
or not with condition (1.27). 
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• In this paper, we propose a method for numerical integration for Riccati 
equation (1.22) which maintains condition (1.27) for each time step At > . We 
detail in second paragraph the simple case of scalar Riccati equation and prove 
the convergence for this particular case ; under some constraints on parameters, 
the scheme is monotonous and remains at the order one of precision, as suggested 
by results of Dieci and Eirola [DE96]. We present the homographic scheme in 
the general case in section 3 and establish its principal property : for each 
time step and without explicit constraint on the time step At, the numerical 
scheme defines a symmetric positive definite matrix. We propose and present 
four numerical test cases in section 4. 

2) Scalar Riccati equation. 

• When the unknown is a scalar variable, we write Riccati equation under the 
following form : 

(2.1) ^ + kx"^ - 2ax - q = 0, 
with 

(2.2) k>0, q >0, 
and an initial condition : 

(2.3) x{0) = d, d > 0, + g2 > 0. 

We approach the ordinary differential equation (2.1) with a finite difference 
scheme of the type proposed by Baraille [Ba91] for hypersonic chemical cinetics 
and independently with the "family method" proposed by Cariolle [Ca79] and 
studied by Miellou [Mi84]. We suppose that time step At is given strictly 
positive. The idea that we have proposed in [Du93], [DS95] and [DS2k] is to 
write the approximation a^j+i at time step {j + l)At as a rational fraction of 
Xj with positive coefficients. We decompose first the real number a into positive 
and negative parts : 

(2.4) a = a+ - a~ , a+ > , a~ > , a+ a~ = , 
and factorise the product x"^ into the very simple form : 

(2.5) (x^).^,,. = x.x.^, . 

^ ^ ^ + 1/2 3 3 + 1 

Definition 1. Numerical scheme in the scalar case. 

For resolution of the scalar differential equation (2.1), we define our numerical 
scheme by the following relation : 

(2.6) — — r + kx. X. , . — 2a^ x. + 2a~ x. , . — q = . 
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• The scheme (2.6) is imphcit because some equation has to be solved in order 

to compute Xjj^i whereas Xj is supposed to be given. In the case of our scheme 
this equation is linear and the solution Xj^i is directly obtained with scalar 
scheme (2.6) by the homographic relation : 

, , (l + 2a+ At)xi -\- qAt 
(2 7) X • 1 = ^ / J ' ^ 

^ ' ^ kAtXj + (1 + 2a- At) 



Proposition 2. Algebraic properties of the scalar homographic scheme. 

Let {xj)j^-^ be the sequence defined by initial condition : 

(2.8) xo = x{0) = d 

and recurrence relation (2.7). Then sequence ixj)^^-^ is globally defined and 
remains positive for each time step. 

(2.9) Xj > 0, Vj € IN, VAt > 0. 
• If At > is chosen such that : 

(2.10) 1 + 2|a|At - kqAf ^ 0, 

then {xj)j^-^ converges towards the positive solution x* of the "algebraic Ric- 
cati equation" : 

(2.11) kx^ - 2ax - q = 

and explicitly computed according to the relation 

1 



(2.12) X* = -(a + + kq ) . 

• In the exceptional case where At > is chosen such that (2.10) is not 

. n 1 1 1 / N . 1 1 1 + 2a+ At , 
satisfied, then the sequence [Xj)^^^ is equal to the constant ^-^^ tor 

j > 1 and the scheme (2.7) cannot be used for the approximation of Riccati 
equation (2.1). 

Proof of proposition 2. 

• The proof of the relation (2.9) is a consequence of the fact that the recurrence 
relation (2.7) defines an homographic function f : 

(2.13) Xj^i = f{xj) 

'yx -\- 

with positive coefficients a , /3 , 7 , : 
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(2.15) a = 1 + 2a+At, ^ = qAt, 7 = k At , ^ = 1 + 2a~ At . 

Because Xq — d > , it is then clear that Xj > for each j > and relation 
(2.9) is established. The homographic function /(•) is a constant scalar equal 

a ,5 _ _ _ 
to — = — when the determinant : 

7 

(2.16) A = a6 - 

is equal to zero. It is the case when condition (2.10) is not realised. When (2.10) 
holds, we have simply : 

(2 17) IM. = - ^7 

^ ■ ^ fix) i-fx + d){ax + (3) 

and / is a bounded monotonic function on the interval ]0,cx)[. Let x* be the 
positive solution of equation : 

(2.18) fix) = X. 

It is elementary to observe that x* is given by the relation (2.12). Let x- be 
defined by : 

(2.19) ._ = -^=..-2^^, 

the other root of equation (2.18). We set : 

(2.20) u, = ^-^ . 
Then we have : 

ax* -\- P axj -\- P 
^ ^x* -\- 5 '^Xj + 6 X- + 6 Xj — X* 

\ ■ ) Uj+i - q;^^. _^ p ^ ax- + (3 ~ 'yx* + S X- - Xj ' 

'jXj -\- 6 jx- -\- 6 

if (2.10) holds. Replacing X- by its value given from (2.19) and using (2.21), 
we obtain : 

6x*-p 

(2.22) Uj+i = Uj . 

^ ' ax* + /3 ^ 

Hi -1-1 

The sequence [uj)^^^ is geometric and the ratio — — has always an absolute 

value less than 1. Effectively we have on one hand : 

-[ax* + 15) < Sx* -P 

because a -\- S = 2 (1 + |a|At) is strictly positive and x* is positive and 
on the other hand : 
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6x* -13 < ax* + 13 
because {a - 6)x* + 2/3 = 2At {ax* -\- q) = At {k{x*Y + q) which is a 
positive number. The absolute value of is not exactly equal to 1 because 



a:* > and ax* + q = x* {kx* — a) = x* ■sj a^ + kq > according to 
(2.10). Then uj is converging to zero and Xj toward x* that completes the 
proof. □ 



Theorem 1. Convergence of the scalar scheme. 

We suppose that the data k,a,q of Riccati equation satisfy (2.2) and (2.10) and 
that the datum d associated with the initial condition (2.3) is relatively closed 
to X* , ie : 

(2.23) --^ + r]<d-x*<C, 

k T 

where C is some given strictly positive constant (C > 0) , x* calculated accord- 
ing to relation (2.12) is the limit in time of the Riccati equation, r is defined 
from data k, a, q according to : 

(2.24) r = , ^ 

2 + 

and rj is some constant chosen such that 

(2.25) < 7/ < . 

k T 

• We denote by x{t;d) the solution of differential equation (2.1) with initial 
condition (2.3). Let {Xj{At ; d^))(^j^-^^ be the solution of the numerical scheme 

defined at the relation (2.7) and let d^ be the initial condition : 

(2.26) xo{At- d^) = d^. 

We suppose that the numerical initial condition d^ > satisfies a condition 
analogous to (2.23) : 

(2.27) --^ +ri < d^- X* < C, 

K T 

with C and rj > equal to the constant introduced in (2.23) and satisfying 
(2.25) . 

• Then the approximated value {xj{At; d^))j^-^ is arbitrarily closed to the 

exact value x{jAt] d) for each j as At — y and d^ — > d. More precisely, 
if a ^ we have the following estimate for the error at time equal to jAt : 

(2.28) \x{jAt; d) - Xj{At; d^)\< A{At -\- \d - d^\), Vj G IN, < At < 5 
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with constants ^ > 0, > , depending on data k, a, g, 77 but independent on 
time step At > and iteration j . 

• If a = , the scheme is second order accurate in the following sense : 

(2.29) \x{jAt]d)-x.{At]d^)\<A{At'^+\d-d^\), Vj€]N,0<At<5 
with constants A et B independent on time step At and iteration j . 

Remark. 

• A direct application of the Lax theorem (see e.g. [La74]) for numerical 
scheme associated to ordinary differential equations is not straightforward be- 
cause both Riccati equation and the numerical scheme are nonlinear. Our proof 
is based on the fact that this problem is algebrically relatively simple. 

Lemma 1. 

Let d^ be some discrete initial condition and x{t d^) be the exact solution of 
Riccati differential equation (2.1) associated to initial condition (2.26). We set 

(2.30) y{t; d^) = x{t- d^) - a;*. 
Then we have 

(d^ - a:*) e-V^ 

(2.31) y{t]d.) = 

1 + kT{d^ - x*){l - e-^n) 

Proof of Lemma 1. 

• The real function IR 9 t 1 — )■ y(t; d^) E JR introduced in relation (2.30) 

satisfies the equation — ^ + ky^ + —y = , that can be rewritten under 

at T 

the form : 

- - ^-^dy) = -dt. 



ky'^ + - ^ y 1 + ^^y 

T 

After integration, we have : 

^a) ^ 1/(0; c^a) -t/r 

1 + kry{t,dj ~ 1 + kry{0,dj ^ 
giving simply : 

« - X*) e-V^ 



1 + kT{d^ - X*) (1 - e-V'^) 
i.e. relation (2.31). Then Lemma 1 is established. □ 
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Lemma 2. 

Let be some discrete initial condition associated to (2.26) and e^. the er- 
ror between the exact solution x{jAt; d^) of the differential equation and the 
solution of the numerical scheme {At ; d^) : 

(2.32) e. = x{jAt;dJ - Xj{At;dJ. 

Let X* and x^ defined in (2.12) and (2.19) be the two fixed points of the 
homographic function /(•) introduced in (2.14) and (2.15). We set 

(2.33) Mf) = ""W"- ^ , (> -1 

and we introduce on one hand the function t i — y z{t) relative to the exact 
solution : 



(2.34) z{t) = 



x{t ; d^) — X* 



tic — tic ^ ^ ^/\ ^ 

and on the other hand a new sequence uj by the same algebraic relation 

X. (At : d.) — X* 

(2-35) Uj = — . 

^ ^ ' X--x.{At;dJ 

Then we have the following estimate : 
(2.36) \e.\ < \h'{z{0))\ \uj - z{jAt) \ . 



Proof of Lemma 2. 

• We have constructed function h{») in order to have h{z{6)) = z{h{6)^ , 
\/e eJR. Then we have 

x{t; d^) = h[z{t)) for each t > 

Xj{At,d^) = h{uj) for each J > 0, 

with function z{») introduced at relation (2.34) and sequence Uj in (2.35). Then 
can be rewritten with the help of this function /i(«) and we have : 

(2.37) e. = h(z{jAt))-h{uj) < sup \h' (z{^)) \\z{jAt) - Uj\ . 

^ e [zijAt) , uj] 



• We note that due to (2.22), the sequence uj is a geometric converging one, 
then we have : — \uo\ < uj < \uo\ for each j > . Moreover thanks to 
relation (2.31) of Lemma 1, we have the following calculus : 
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y(t : a 

z{t) 



X- -X* -y{t; cL ) 



A' 

'A 



1 + kT{d^ - x*){l - e~^/^) 



{d^ - X*) e-^/^ 



1 + kT{d^ - x*){l - e-^/^) 

krjd^ - X*) 

~ 1 + kT{d^ - X*) ^ 

thus satisfy clearly : — |2;(0)| < z{jAt) < \z{0)\ , Vj > 0. We remark 
also that : 

krld. — X*) 

-1 < z(0) = no = -T^T^. 

1 + kr [d^ — X*) 

• Morever, is a decreasing function for ^ > — 1. Then by the two points 

finite difference Taylor formula and the above statements, we deduce from (2.37) 
the estimate (2.36) and Lemma 2 is established. □ 



Lemma 3. 

Let r be defined in (2.24). We introduce a and (3 according to the following 
relations 

(2.38) a = ^ - \a\ 



2t 
1 



(2.39) /3 = -+\a\ 

and when At > we define function ^p{At) according to 

(2.40) ,(At) = l + ^log(i^). 

Then with notations introduced in Lemma 1 and the following one 

(2.41) % = ^, 

we have 

(2.42) Uj - z{jAt) = -z{0) e 1 - exp 
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Proof of Lemma 3. 

• We study the difference that majorate the right hand side of (2.37). From 
relation (2.21), we know that sequence Uj is geometric and more precisely : 



Uj = 



jx* + S 

with 7 and 6 computed in (2.15) j = k At , (5 = 1 + 2a~ At . We have 
from equation (2.11) and defintion (2.24) of variable r : 



+ S 



kAt ^ (o - + 1 + 2a" At 



^x* + d 



kAt I + ^ ) + 1 + 2a" At 
At (^a+ - c 



2r 

- —] + 1 + 2a~At 



2tJ 



At (a+ 



a' 



+ ^] + 1 + 2a~ At 



1 



2tJ 



1 + At{\a\ --) 

1 + At(|a| +^) 
1 - aAt 



1 + l3At 
We can now write : 



z{jAt) 



= -z{0) e 



-z(0) e" 



-jAt/r 



-jAt/r 



1 - 



1 - aAt\i 
1 + /3At 



/I - aAt 

exp I ■ + ? log — — 

^\ 3 ^Vl + 13 At 



exp {^j 



T 



1 + At'°S 



aAt 



1 + ^At 




and relation (2.42) is proven. Lemma 3 is established. 



□ 



Lemma 4. 

Let (^(•) be defined by relation (2.40). We suppose that time step At satisfies 
the condition 

(2.43) < At < ^. 

Then we have 
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(2.44) 



ifiiAt) 



T 



Proof of Lemma 4. 

• We have the following elementary calculus 



1+^ 

and we deduce after integration 



if i^i< 1 



< 



iog(i + - - ^ + C 



< -— + — < — 
45-2 



if 1^1 < 



Then we have 
(2.45) 



log(l + - - ^ + C 



< -1^1 
- 2 



if i«i < 2 



and we use this estimation with ^ = —a At and ^ = (3 At. We remark that 

I a I < — then (3 < — and we deduce that when condition (2.43) is satisfied, 
2 r r 

we have 

At < - < — < — , then aAt < - and /3 At < - . 
- 2 - 2/3 ~ 2a ~ 2 ^ "2 

We deduce from (2.45) 

1 - aAt 



log 



1 + ^At 



{-OL - l3)At + 



13' 



< -(a^ + l3'^)At' 



< 



• As a consequence of the previous inequality, multiplying the above expression 

by we obtain : 
At 



<f{At) 



-a - P)At + ^^—^At" + \{-a^ - f)Af 



< 



< -(«^ + /3^)rAt^ 
1 1 

and due to the fact that i—a — /3)r = —1 we have since r < — < — 

p a 



<p{At) - 



and relation (2.44) is proven. Lemma 4 is established. 



□ 
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Lemma 5. 

We suppose that a 7^ 0. If time step At satisfy 
(2.46) < At < inf |a| r^), 

we have 



(2.47) I <f{At) 



7 

\a\ At \ < — \a\ At 
II I - 12 ' ' 



Proof of Lemma 5. 



If a ^ 0, we observe that : ,5^ — o;^ = 



2\a\ 



T 



> , and inequahty (2.46) 



suppose that At has been chosen such that At < | a | . Due to previous 
computations, we have the following set of estimations : 



a 



r At - ip{At) 



< {^(a^ + /3')At + i(a4 + /34)At2) rAt 



At / 1 3a 
I + 



< 



3 V4r3 r 
At At^ \ 
3^3 ^ 2^) 



2^ At2 

+ -— 



) 

r At 



1 



r At 



< 



< 



< 



At^ /I 1\ 
V2 



because |a| < — due to (2.24) 



7 At^ 



12 



and relation (2.47) is proven. 



□ 



Lemma 6. 

We suppose that a = 0. If time step At satisfy condition (2.43), we have 



(2.48) 



(^(At) + 



12 r2 

Proof of Lemma 6. 



At^ 



< 



32 



At^ 



If a = 0, then (3 = a = — and we have simply 



-(a^+^^)rAt^ 



2 / 1 \3 .2 
- — rAt^ 



3 V2t 



1 



12 



At^ 
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-(a^ + ^^)rAt3 = (—YrAt^ = At^ < At^ due to the 

2 ^ ^ ^ \2tJ 16r3 - 32r2 

relation (2.43). Then the relation (2.48) is a direct of previous estimation (2.44) 
established in Lemma 4. That completes the proof of Lemma 6. □ 



Lemma 7. 

We suppose that function t i — > z{t) is defined at relation (2.34) and that the 
numerical initial condition satisfies hypothesis (2.27). We denote by h{») 
the expression introduced in (2.33). Then we have 

(2.49) I z(0) I < - mm(c , -T]] 

7] \ kr / 

(2.50) \h'{z{0))\ < kr (min (c , —-ri^^ . 
Proof of Lemma 7. 

• We have, due to relation (2.27) : —C < x* — d. < - — — ij . 

^ kr 



—C — ^ X* ~ TZ ~ = X- ~ ^ —"U in consequence 



Then following (2.19) and (2.24) we deduce 

kr kr 

1 1 
2.51 I — < -. 

\x- — d/^] 7] 

From relation (2.34) we have 2;(0) = — and inequality (2.49) is a 

X- — ga 

direct consequence of (2.51) and of hypothesis (2.27). We derive now the ex- 
pression (2.33) relatively to variable ^ and we have easily : 

due to the above expression of 2;(0) and relation (2.51). The estimation (2.50) 
is established and Lemma 7 is proven. □ 

Proof of Theorem 1. 

• We cut the expression inside the absolute value of left hand side of (2.28) 
into two parts. The first one is the error for the continuous solution of the 
ordinary differential equation when changing initial data and the second one is 
to the error t. between the solution of the ordinary differential equation and 

the discrete solution given by the scheme for the same initial condition : 
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\x{jAt;d) - Xj{At; dj\ < \x{jAt- d) - x{jAt; dj \ + \e.\ 
and due to definition (2.30), 

(2.52) I x{jAt ; d) - Xj{At ; d^)\ < \ y{jAt ; d) - y{jAt ; d^) | + | e. | . 

• We first study the term | y{jAt; d) — y{jAt; d^) \ in riglit hand side of 

(2.52) . We first minor ate the absolute value of the denominator of the right 
hand side of (2.31). Under the hypothesis (2.27), we have : 

1 + kT{d^ - x*){l - e~^/^) > 1 if - ar* > 0, 

and in the other case when d. — x* < then \d. — x*\ < - — — rj due to 

^ ^ kr 

(2.27) and we have in consequence 

1 + A;r(d^ -a;*)(l -e~^/^) > 1 - {l-kr 7]) = kr 7] > d^ - x* < . 

Therefore the denominator of right hand side of (2.31) is always strictly positive 
and is in all cases minor at ed by krr] : 

1 + kT{d^ - x*){l - e~^/'^) > krr] Vt>0. 

In consequence of (2.33) and previous algebra, y{t; d) and y{t; d^) are well 
defined for each t such that < t < oo and we have : 

x{jAt; d) - x{jAt- d^) = y{jAt; d) - y{jAt; dj 

{d-d^)e-'^l^ 

{1 + kr{d - x*){l - e-V^)) (l + kT{d^ - x*) {1 - e^V^)) ' 

then : 

(2.53) \xijAt,d) - x{jAt,d^)\ < \d-d^\ 

and the estimation of the first term introduced in (2.52) is then controlled . 

• We now study the error e^. . Due to lemmas 2 and 3 and in particular 

relations (2.36) and (2.42), we have 

(2.54) |e.| < 1^(0)1 \h'(ziO))\ e'^^J | 1 - exp(^,- ^(At) ) . 

j At 

(i) If a = 0, then (p{At) is negative due to relation (2.48) and 6j = 

r 

remains positive. We deduce that in this case 
1 - exp(^ej(p{At)^ < 1 - exp(^ej (p{At)^ < \ej(p{At)\= Oj \(piAt)\ 
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and in consequence of (2.54) and (2.48), 



t2 



< 1^(0)1 \h'(z{0})\ (e-«J e/) (- + -) 

< ^ |. (0)11;.' (40)) I ;i^Pp ^ 

and relation (2.29) is established with 



(/cr?]) ' 96e r?7 ( fcr ^)) 



(2.55) A = inf 

(2.56) B = 

(ii) If a ^ 0, then ^{At) is positive due to relation (2.47) if At satisfies 
condition (2.46). We suppose moreover that time step satisfies also the condition 

(2.57) At < ^ 



19 |a| 

and due to (2.47), we have 

(2.58) ^{At) < (1 + 1) |a| At < i. 

In order to major ate the expression e~^j 1 — ex.p(^Oj (p(At)^ 

between two cases. On one hand, when Oj ^{At) < 1 , we have 
the exponential function over the interval [0, 1] : 

< e^j'^(^*) - 1 < {e-l)ej^{At) and we deduce 



we distinguish 
3y convexity of 



/ \ Q 19 e — 1 
1 - exp(^j(^(At)) < {e-l)[eje~^j]if{At) < — \a\ At . 

The previous inequality and estimation (2.54) show that under hypotheses (2.46) 
and (2.57) concerning the time step, we have 

1 3 

(min (C , 77)) At due to (2.49) and (2.50) 



19 e — 1 a kr 
< — 



12 e 7] 
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19 e - 1 k 

< 

24: e 7] 



- (min [c , —-rj^j At due to (2.24) 



and 

19 1 A; / / 1 \\3 
(2.59) |e^. I < —- -[min^C , — -r]jj At when Oj ^{At) < 1 

On the other hand when 9j (p{At) > 1 , we have : 
1 - exp(ej^{At)j = 
= e^iMAt)-l/2) 

< 
< 



1 - exp -Oj (f{At) 



-Oj/2 ^ _ exp(-6'j(^(At)) due to (2.58) 

3/ Oj (p{At) because 6j and (p{At) are both positive 



< 2 sup (Oe ^) (p{At) < 
0>0^ ^ 



19 2 
12 e 



\a\ At 



thanks to relation (2.58). Following inequality (2.54), we obtain in this second 
case 

7]]] — \a\ At 



T] 



^min 



kr 



6 e 



19 A- / / 1 \\3 

(2.60) |e^. |< j^-[mm{C,—-rijj At when ejip{At)>l 

and relation (2.29) is proved for this case with 

/ / 1 \2 IQ k 

(2.61) A = inf ^ ^ ^ 



(2.62) B = ini ( ^ , \a\ t\ ^ ^ 



19 |a| 



which ends the proof of Theorem 1. 



□ 



3) Matrix Riccati equation. 

• In order to define a numerical scheme to solve the Riccati differential equa- 
tion (1.22) with initial condition (1.19) we first introduce a real number fi, which 

is chosen positive and such that the real matrix [/x/ — {A + A^)] is definite 
positive : 

(3.1) > , ^{fix , x) — {Ax , x) > , \/x ^ . 
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Then we define a definite positive matrix M wich depends on strictly positive 
scalar ji and matrix A : 

(3.2) M = ^iil - A. 

The numerical scheme is then defined in analogy with relation (2.6). We have 
the following decomposition : 

(3.3) A = A+ - A- , 
with 

(3.4) = /// , A~ = M, > 0, M positive definite . 

Taking as an explicit part the positive contribution A'^ of the decomposition 
(3.3) of matrix A and in the implicit part the negative contribution A~ = M 
of the decomposition (3.3), we get the following harmonic scheme : 

(3.5) I Xt^^^+' ~ + li^jKXj+i + Xj+,KXj) + 

The numerical solution Xjj^i given by the scheme at time step {j + 1) At is then 
defined as a solution of Lyapunov matrix equation with matrix X as unknown : 

(3.6) S^jX + XSj = Yj. 
with : 

(3.7) Sj = h + ^KXj + AtM, 
and : 

(3.8) Yj = Xj -\- fiAtXj + AtQ. 

We note that Sj is a (non necessarily symmetric) positive matrix and that Yj is 
a symmetric definite positive matrix if it is the case for Xj. In all cases, matrix 
Yj is a symmetric positive matrix. 

Definition 2. Symmetric matrices. 

Let n be an integer greater or equal to 1. We define by <S„(IR), (respectively 
4S+(IR) , <S+*(IR) ) the linear space (respectively the closed cone, the open cone) 
of symmetric-matrices (respectively symmetric positive and symmetric definite 
positive matrices) ; we have 

(3.9) {x, Sy) = {Sx,y), Va:,2/€lR", VS'€5„(IR), 

(3.10) {x,Sx) > 0, yxeW, V5'€5+(1R), 

(3.11) {x,Sx)>0, Wxen'^^x^O, VS'e<S+*(IR). 
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The following inclusions «S+*(IR) C <S+(IR) C <Sn(IR) are natural. 

Proposition 3. Property of the Lyapunov equation. 

• Let 5 be a matrix which is not necessary symmetric, such that the associated 
quadratic form: IR"" 3 x i — )■ {x^ Sx) € IR, is strictly positive : 

(3.12) S + g5+*(IR). 

Then the application ip^ defined by : 

(3.13) 5n(IR) 3 X ^ ipg{X) = S^X + XS e 5n(IR), 

is a one to one bijective application on the space <S„(IR) of real symmetric 
matrices of order n. 

• Morever, if matrix (Pg{X) is positive (respectively definite positive) then 
the matrix X is also positive (respectively definite positive) : 

i^fgiX) G <S+(IR) =^ X e S+m) and (^^(X) G <S+*(IR) ^Xe <S+*(IR)) . 
Proof of proposition 3. 

• We first observe that (p^ is a linear map. In the case where 5* is a symmetric 
matrix, we can immediatly deduce from (3.12) that S is definite positive. Let 
X be a matrix such that : 

(3.14) ^^{X) = 0. 

We prove that X — 0, ie that kenp^ = {0}. Let x be an eigenvector of 
matrix S associated to eigenvalue A : 

(3.15) Sx = Xx, x^O. 

We deduce from relations (3.14) and (3.15) and the fact that matrix S is sup- 
posed to be symmetric the equality 

(3.16) S{Xx) = -X{Xx). 

According to the previous hypothesis the negative scalar —A cannot be an eigen- 
value of S , then the vector [X x) must be equal to zero. This is true for each 
eigenvalue of matrix S, wich prove the property in this case, because (p^ is also 
an endomorphism from tS„ (IR) to Sn (H) . 

• In the case where S is not symmetric we suppose first that S is composed 
into blocks of Jordan type. We first study the case where 5" = A is composed 
of exactly one Jordan block associated to eigenvalue A : 
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/A 





(3.17) A 



1 
A 




1 









••• A 1 
\0 ••• A/ 

and according to the hypothesis (3.12), we have : 

(3.18) Re A > 0. 

Let Xij be an element of the matrix X , we have from the following relation 

ifi > 2 andj > 2 



{A^X + XA) 



2AXij -\- Xi- 



(3.19) 
and : 
(3.20) 

(3.21) 

Because A 7^ 0, Xi^i = from (3.20), then Xij 
By induction, Xij — taking into account (3.13). 



{A^X + XA)i,i = 2AXi,i 
{A^X + XA)i,, = 2AXi,, + 



J >2. 

if i > 2 from (3.21). 



• When matrix S is composed by a family of Jordan blocks of the previous 
type, i.e. S = A = (diagAj) where 1 < j < p and Aj is a Jordan block of 
the type (3.17) and of order rij 

••• \ 
••• 



/Ai 




(3.22) A = 





A2 





V 



A 







p-i 






A J 



we decompose the matrix X into blocks Xij of order rii x rij 



( 



(3.23) X = 



^1,1 

A^2,l 



Xi^2 
A^2,2 



-^1,3 



X 



2,3 



X 

X 



2,P 



\ 



X 



\ X 



p-1,1 



X 



X. 



P,2 



p-l,p-l 
Xp^p-i 



X 



p-l,p 



X. 



I 



Then the block number (i,ji) of the expression S'^ X-\-XS is equal to A^X^ j + 
Xi j Aj and for i ^ j we have to prove that the nondiagonal matrix X 



identically null. 



• We establish that if A is a Jordan block of order n of the type (3.17) 
satisfying inequality (3.18), if M is a second Jordan block of order m of the 
previous type. 



22 



HOMOGRAPHIC SCHEME FOR RiCCATI EQUATION 



f II 1 ••• ON^ 
1 ••• 



(3.24) M = 



••• fi 1 
\0 ••• fi/ 

such that an inequahty analogous to (3.18) is valid : 

(3.25) Refi > 0. 

and if X is a real matrix of order n x m chosen such that 

(3.26) A^X + XM = 0, 

then matrix X is identically equal to zero. As in relations (3.19)-(3.21), we 
have : 

(3.27) {A^X + XM)i,i = AXi,i + 

(3.28) {A^X + XM)ij = XXi^j + /iXij + j > 2 

(3.29) {A^X + XM),,i = AX,,i + + z > 2 

(3.30) (A^X + XM)ij = XXij + //Xij + A,_i,, + A,j_i ^, j > 2. 

Then due to (3.18), (3.25) and (3.27) we have Ai i = 0. By induction on j and 
according to relation (3.28) we have Xij = 0. By induction on i and due to 
relation (3.29) we have analogously A^^i = 0. Finally relation (3.30) prove by 
induction that Xij = when i and j > 2, and matrix X is identically null. 

• When matrix S = A = (diagAj) is composed by a family of Jordan blocks 
of the previous type, the nondiagonal blocks Xij of decomposition (3.23) are 
null due to the previous point. Moreover, the diagonal matrices Xi^i are also 
identically null due to the property established at relations (3.17) to (3.21). 
Then the proposition is established when matrix 5* = A is composed by Jordan 
blocks as in relation (3.22). 

• In the general case where real matrix S satisfy relation (3.14) there exists 
always a complex matrix Q such that : 

(3.31) S = Q-^AQ, 

where matrix A has a bloc Jordan form given e.g. by the right hand side of the 
relation (3.22). We have also the following elementary calculus : 

iPg{X) = S^X + AS' 

= A^ X + XQ-^ AQ 

= A^ {Q-^XQ-^)Q + {Q-^XQ-^)AQ 

= (A* Y + YA)Q 
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with Y = Q~^XQ~^ . Because the matrix Q is inversible and ^g{X) is equal 
to zero, the matrix (A^ Y + YA) is equal to zero. Matrix A is a diagonal bloc 
Jordan form, and from the previous points we deduce that Y is equal to zero, 
then X = and the proof of first assertion of proposition 3 is established in 
the general case. 

• We suppose now that the matrix ip^ {X) is symmetric and positive, that is : 
^g{X) G tS+(IR). Let X be an eigenvector of matrix X associated to the real 
eigenvalue A : 

(3.32) Xx = Xx, X ^ 0. 

From the definition of application (f^ , we have the following relation : 

(3.33) {x,cp^{X)x) = 2X{x,Sx) 

and from hypothesis on (f^ {X) the left hand side of relation (3.33) is positive. 
The expression {x, Sx) is also strictly positive since relation (3.8) holds. We 
deduce that the number A is positive because X G tSri(]R) has an orthogonal 
decomposition in eigenvector spaces. 

• If matrix ^g{X) is symmetric and positive definite {ipg{X) G tS^*(IR)), 
we introduce eigenvector x of matrix X as previously (relation (3.32)). Then 
relation (3.33) remains true and the left hand side of this relation is strictly 
positive. Then the eigenvalue A of matrix X is strictly positive and proposition 
3 is proven. □ 



• The numerical scheme has been writen as an equation of unknown X = 

Xj^i that takes the form : 

(3.34) ^^^{X) = Yj. 

with (p^ given by a relation of the type (3.13) with the help of matrix Sj 
defined in (3.7) and Yj in (3.8). Then we have the following proposition. 

Proposition 4. 

Homographic scheme computes a definite positive matrix. 

• The matrix Xj defined by numerical scheme (3.5) with the initial condition 

(3.35) Xo = 

is positive for each time step At > : 
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(3.36) XjeS+{M), Vj>0. 

• If there exists some integer m such that Xm belongs to the open cone 
<S^*(IR) , then matrix X^+j belongs to the open cone <S^*(IR) for each j : 

(3.37) (3m G IN, G <S+*(IR)) =^ (vj > 0, X^+j G 5+*(IR) ) . 
Proof of proposition 4. 

• First we have Yq = AtQ and .So = ^/ + AtM, then Xi is a symmetric 
positive matrix {X G tS+(IR)) according to proposition 3 since matrix 5*0 is 
symmetric positive and M has been chosen such that 

(3.38) M + M * is positive definite. 

The end of the first point follows by induction. 

• If real symmetric positive definite matrix Xjj^i is given, relation (3.8) clearly 
indicates that matrix Yj is symmetric positive definite and matrix Sj intro- 
duced at relation (3.7) has a symmetric part which is positive definite if we 
verify that the following matrix 

(3.39) KXj + XjK 

is positive. But this property is a consequence of the following. On one hand, 
matrix K is positive definite and map X i — y K X + X K transforms the 
closed cone of positive symmetric matrices onto himself. On the other hand, 
matrix Xj is symmetric positive definite then K Xj + XjK G tS+*(IR) due 
to proposition 3. Proposition 4 is established. □ 

• We have defined a numerical scheme for solving in a approximate way the 
Riccati equation (1.20) with the help of relation (3.5) and the initial condition 

(3.40) Xo = 

naturally associated with initial condition (1.19). Recall that time step At is 
not limited by any stability condition : matrix Xj is allways symmetric posi- 
tive and positive definite if matrix Q is symmetric positive definite. Moreover 
the equation (3.5) that allows the computation of Xjj^i from data is a lin- 
ear equation whose unknown is a symmetric matrix. But Xj^i is a nonlinear 
(homographic !) function of previous iteration matrix Xj. 

• We now study the convergence of the iterate matrix Xj as long as discrete 
time j At tends to infinity. We know from Proposition 1 that the solution of the 
differential Riccati equation tends to the solution of stationary Riccati equation : 

(3.41) XooKXoo - {A^Xoo + XooA) - g = . 



25 



FRANgois Dubois and Abdelkader Saidi 

We first study the monotonicity of our numerical scheme. Recall first that if A 
and B are two real symmetric matrices, the condition 

(3.42) A < B 

and respectively the condition 

(3.43) A < B 

means that matrix B — A is positive (^B — A E «S+(IR)) 

(3.44) {x,{B-A)x) > 0, Vx-elR^, 

and respectively that matrix B — A is positive definite (^B — A E «S^*(IR)) : 

(3.45) {x, {B-A)x) > , Va:€lR". 

Proposition 5. Monotonicity. 

• Under the two conditions : 

(3.46) Q is a definite positive symmetric matrix 

(3.47) i (KXoo + X^K) < (/^ + ^) ^ , 

the scheme (3.5) is monotone and we have more precisely : 

(3.48) (o < Xj < Xoo) ^ (O < Xj < Xj+i < Xoo). 

Proof of proposition 5. 

• We know from Lewis [Le86] that for symmetric definite positive matrix K 
and symmetric positive matrix Q, the algebraic Riccati equation (3.41) has 
a unique symmetric positive solution X^o- Moreover, matrix Xoo is positive 
definite if matrix Q is positive definite. 

• We first establish that matrix O = Xq^ — Xj^i is positive if the left hand 
side of implication (3.48) is satisfied. We substract relation (3.41) from numerical 
scheme (3.5), observe that 

Xj K Xj^i — Xqq K XoQ = Xj K [Xj^i — Xoo) + {Xj — Xqq) K Xqq 

= {Xj-Xoo)KX^ - XjKe 

and that 

-^3+1 -^3 ~ Xqq K Xqq = {Xj^i — X^) K Xj + Xoo K [Xj ~ Xqq) 

= X^K{Xj - Xoo) - QKXj 
then we obtain the following equation for G : 
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(3.49) < 



+ i {XjKe + SKXj) + (M^e + eM) = 

(Xoo - Xj) + M (Xoo - X,) + i [(X,- - Xoo) i^Xc 
+ Xoo K {Xj - Xoo)] = {Xao - Xj) 



with 
Si = 



^/ + -KX^ + M 
2 At 2 ^ 



The matrix Ei admits a positive definite symmetric part Ei + because it is 
the case for matrix /. Moreover, since matrix Xoo — Xj is symmetric positive 

by hypothesis, it is sufficient to estabhsh that the symmetric matrix E2 + Tj\ 
is positive definite and to apply the proposition 3. This last property is exactly 
expressed by hypothesis (3.47) and the first point is proven. 



We consider now the matrix 5*00 defined by : 



(3.50) 



Sr 



= -{XooK + KX^) - (A* + A) 



The matrix 6*00 is symmetric and we have the following calculus 
1 

00 ' 00 



5, 



00 



(Xoo K Xoo) + X-' (Xoo K Xoo) 



X-l(XooA) + (A^Xoo)X, 



-1 

00 



(3.51) 
with 



Soo = (f^ (Xj-) 



.t 



00 



A^ X. 



00 



Then due to relation ( 3.41) we 



have E3 + S3 = Q > due to hypothesis (3.46). Then the propostion 
3 joined with (3.51) and the fact that matrix X^^ is symmetric positive definite 
establish that matrix ^oo is symmetric definite positive. 



• We establish that under the same hypothesis (3.48), the matrix Z = Xj^i — 
Xj is positive. We start from the numerical scheme (3.5) and replace matrix Q 
by its value obtained from relation (3.41). It comes : 

^^^Z) ^ ^^Z + \{XjKZ + ZKXj) + {M^Z + ZM) = 
= fiXj + g - XjKXj - {M^Xj + XjM) 
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= Xj + XjA + Q - XjKXj 

= -'^oo K Xao — Xj K Xj — [A^ [Xoo — Xj ) + [Xoo 

K {Xoo — Xj) + {Xqq — Xj) K Xoo + Xj K {Xoo — Xj) + 



X,)A] 



1 
2 



+ (Xoo - Xj)KXj - [A^ (Xoo - Xj) + (Xoo 



= if^ [Xoo - Xj) with S4 



]^{KX^ + KX,) 



Xj)A] 
- A. 



The matrix 



S4 + S5 = - [XjK + KXj) + 5c 



is positive definite due 



to the previous point ; in consequence matrix ip^^ (-^00 — Xj) is symmetric 

positive. The end of the proof is a consequence of propositon 3 and of the fact 

1 1 

that the matrix Si = -I+-KXj+M has clearly a symmetric part 



El + which is positive definite. 



□ 



Proposition 6. Convergence when discrete time tends to infinity. 

We suppose that the data A^ Q of Riccati equation (1.20) and parameters 
and At of harmonic scheme (3.5) satisfy the conditions (3.1), (3.46) and (3.47). 
Let Xj and X^o be the solution of scheme (3.5) at discrete time j At and 
the symmetric definite positive matrix solution of the so-called algebraic Riccati 
equation (3.41). Then Xj tends to X^q when j tends to infinity : 

(3.52) Xj X^. 
Proof of proposition 6. 

• Let £k be the Grassmannian manifold composed by all the linear subspaces 
of space IR" with dimension exactly equal to k : 

(3.53) £k = {W,W subspace of IR", dimW = k} . 

Then we have the classical characterization of the k° eigenvalue of symmetric 
matrix A with the so-called inf-sup condition (see e.g. Lascaux and Theodor 
[LT86]) : 

( Av 

(3.54) Lbk = inf sup —, — . 

We£k veW (^'^) 

Let Xj be the k° eigenvalue of matrix Xj, the k° eigenvalue fik (yU-i < 
/i2 < • • • < fJ'n) of matrix Xq^ and W a fixed subspace of IR"^ of dimension 
exactly equal to k. We deduce from inequality (3.48) of Proposition 5 : 
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{Xj V,V) {Xj+i V,V) (Xoo V , v) 



< 



< 



Then when subspace W is arbitrarily given in set we have : 



sup 



iy, v) 



< sup 

veW 



iy, v) 



< sup 

veW 



(Xqo V , V) 

{v, v) 



and taking the infimum bound of previous hue as subspace W belongs to Grass- 
mannian manifold £k we deduce, thanks to (3.54) 

(3.55) A,^ < A,\i < . 

By monotonicity, eigenvalue A^ is converging towards some scaler fj,^ as j tends 
to infinity : 



A 



.k 



> /I as J 



oo , 



l<k<n 



(3.56) 

• Consider now the unitary eigenvector Vj of matrix Xj associated with 
eigenvalue Xj : 

(3.57) Xj = ^jVj , II II = 1 , l<k<n, j>0. 

Because Xj is a symmetric matrix, the family of eigenvectors {'^j)i^j^<^^ is 
ort honor mal and defines an orthogonal operator pj of space IR"^ defined as 
acting on the canonical basis of space IR" by the conditions 

(3.58) pj.ek = , 1 <k <n, j >0. 

Rotation pj belongs to the compact group 0{n) of orthogonal linear transfor- 
mations of space IR"" . Then after extraction of a convergent subsequence Pj of 

the initial sequence {pj)jyQ we know that there exists an orthogonal mapping 

Poo £ 0{n) such that 

(3.59) p'j — > Poo as j — > +(X) . 

We introduce the family (w^^)i<;[,<^ of vectors in IR"^ by the conditions 



(3.60) 



w' 



= Poo»efc, I <k <n 



It constitutes an orthogonal basis of space IR" and for each integer k{l < k < 



n) , the extracted sequence of vectors Vj^ is converging towards vector w[ 



(3.61) 



V 



Ik 



1 < < n , j — )■ +00 



• We introduce the symmetric positive definite operator Y^o by the conditions 
(3.62) Y^.wi^ = p^w^^ ,l<k<n. 

We study now the convergence of the subsequence of matrices X'j towards Yoo • 
We first remark that the sequence of matrices {Xj^ .^^ is bounded in space 
<S„(IR) : 
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(3.63) II Xj II < , V j € IN 

and we have also the following set of identities : 

{Xj - Yoo)wt = Xjiwi, - v^) - Xto^l + \)v] 

= X,{wl, - v]) - A^(w;^ - v)) + (A^^ - \l,)v). 
So we deduce, due also to (3.63) and (3.57) : 

(3.64) II (X, - Yoo)wl^ II < A^ II wl^ - v] \\ +A^ || w;^ - v] \\ + | 

and according of the convergence (3.61) of subsequences Vj and (3.56) of se- 
quences A^ as index j tends to infinity, we get from estimation (3.64) : 

(3.65) X'j - Yoo ^ 0, j ^ +00 . 

• Due to the definition (3.5) of the numerical scheme, sub-sequence Xj nec- 
essarily converges to the unique positive definite matrix Xoo of the algebraic 
Riccati equation and in consequence we have necessarily 

(3.66) ^00 = ^00- 

We deduce that for any arbitrary subsequence of the family (Xj)^.^^ there exists 

an extracted sub-subsequence converging towards X^o as j tends to infinity. 
Then the entire family (Xj)^.^^ converges towards X^q as j tends to infinity 
and the property is established. □ 



4) Numerical experiments. 
4.1) Square root function. 

• The first example studied is the resolution of the equation : 

(4.1) ^ + x2 - g = 0, x{0) = 

with n = 2, A = 0, K = I and matrix Q equal to 

(-) ^-Ki ")G loo) (-\ ;)■ 

• We have tested our numerical scheme for fixed value At = 1/100 and 
different values of parameter n : n = 0.1, 10-^, 10+^. For small values of 
parameter /x, the behaviour of the scheme does not change between // = 0.1 
and n = 10~^. Figures 1 to 4 show the evolution with time of the eigenvalues 
of matrix Xj and the convergence is achieved to the square root of matrix Q. 
For large value of parameter fi [fi = 10+^), we loose completely consistency of 
the scheme (see figures 5 and 6). 
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Figures 1 and 2. Square root function test. 
Two first eigenvalues of numerical solution (// = 0.1). 





Figures 3 and 4. Square root function test. 
Two first eigenvalues of numerical solution (// = 10~^). 



ftlsft cijgetiyalue 



second dgenvalue 





Figures 5 and 6. Square root function test. 
Two first eigenvalues of numerical solution {/j, — 10+^). 
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4.2) Harmonic oscillator. 

• The second exemple is the classical harmonic oscillator. Dynamical system 
y{t) is governed by the second order differential equation with command v{t) : 

(4.3) ^ + 2S^ + .^yit) = bvit). 

This equation is written as a first order system of differential equations : 

The parameters R, Q, cj, 6 and b of the ordinary differential equation (4.4) and 
the cost function (1.6) are given by : 




Figures 7 and 8. Harmonic oscillator. 
Two first eigenvalues of numerical solution (// = 0.1, a = 0.01, At = 0.01). 

• In this case, we have tested the stability of the scheme for fixed value of 
parameter = 0.1) and different values of time step At. We have chosen 
three sets of parameters : a = At = 1/100 (reference experiment, figures 7 
and 8), a = 10~® , At = 1/100 (very small value for a , figures 9 and 10) and 
a = 1/100 , At = 100 (too large value for time step, figures 11 and 12). Note that 
for the last set of parameters, classical explicit schemes fail to give any answer. 
As in previous test case, we have represented the two eigenvalues of discrete 
matrix solution Xj as time is increasing. On reference experiment (figures 7 
and 8), we have convergence of the solution to the solution of algebraic Riccati 
equation. If control parameter a is chosen too small, the first eigenvalue of 
Riccati matrix oscillates during the first time steps but reach finally the correct 
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values of limit matrix, the solution of algebraic Riccati equation. If time step is 
too large, we still have stability but we loose also monotonicity. Nevertheless, 
convergence is achieved as in previous case. 




Figures 9 and 10. Harmonic oscillator. 
Two first eigenvalues of numerical solution ( = 0.1, a = 10~^, At = 0.01). 



fijttst cijgetiyaltie 



second dgeuvalue 




km™ 




Figures 11 and 12. Harmonic oscillator. 
Two first eigenvalues of numerical solution { fi = 0.1, a = 0.01, At = 100). 



4.3) String of high speed vehicles. 

• This example has been considered by Athans, Levine and Levis [ALL67] in 
modelling position and velocity control for a string of high speed vehicles. Let 
N be some integer and 

(4.6) n = 2N - 1 

be the order of the given matrices An, and Qn- The matrices Ajq 
and Qn admit the following structure : 
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fail 




(4.7) 



An = 



ai2 

^22 



(4.8) 
(4.9) 



Kn 
Qn 



••• 
••• 
\ ••• 

diag (1,0 





^23 










\ 




a 



N-2N-2 





a 



N-2N-1 



a 



'N-lN-l 





1, 









• , 1,0, 1 

diag ( 10 , , 10 , , • • • , 10 , , 10 , , 

The unknown positive definite matrix Xn satisfies the algebraic Riccati equa- 
tion 





-1 
-1/ 

0, 1) 

• • 10 , 



10) 



(4.10) 



Xn Kn Xn — {A% Xn + Xn An ) — Q 



N 







The solution of this equation is detailed on 'figure' 13 with 10 significative dec- 
imals. The first six decimals are absolutly identical to the ones published by 
Laub [La79]. 



columns 1 to 3 : 

+1.3630206938E + 00 
+2.6172154724^ + 00 
-7.0542734123^ - 01 
+9.3685970173^; - 01 
-2.9366643189^; - 01 
+4.7735386064^; - 01 
-1.9737508953^; - 01 
+2.1121165236^ - 01 
-1.6655183115^-01 

columns 4 to 6 : 

+9.3685970173^ - 01 
+1.4752196951E + 00 
+2.1577096313^; + 00 
+8.2576995027^; + 00 
-1.9464979789^; + 00 
+1.7558734280^ + 00 
-6.7071749345E - 01 
+6.6514730720E - 01 
-4.7735386064E - 01 



+2.6172154724E + 00 
+7.5925521955^ + 00 
-1.6803557707^; + 00 
+1.4752196951E + 00 
-4.5950584109^; - 01 
+6.6514730720^; - 01 
-2.6614220828^; - 01 
+2.8065373289^; - 01 
-2.1121165236^-01 

-2.9366643189^ - 01 
-4.5950584109E - 01 
-6.0913599887^; - 01 
-1.9464979789^; + 00 
+1.8056048615^; + 00 
+1.9464979789^ + 00 
-6.0913599887^; - 01 
+4.5950584109^ - 01 
-2.9366643189^; - 01 



-7.0542734123E - 01 
-1.6803557707E + 00 
+1.7747816032^ + 00 
+2.1577096313^; + 00 
-6.0913599887^; - 01 
+6.7071749345^; - 01 
-2.6284317351^; - 01 
+2.6614220828^ - 01 
-1.9737508953E - 01 

+4.7735386064E - 01 
+6.6514730720^ - 01 
+6.7071749345^; - 01 
+1.7558734280^; + 00 
+1.9464979789^; + 00 
+8.2576995027E + 00 
-2.1577096313^; + 00 
+1.4752196951E + 00 
-9.3685970173^; - 01 
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columns 7 to 9 : 

-1.9737508953E' - 01 
-2.6614220828^ - 01 
-2.6284317351E - 01 
-6.7071749345E - 01 
-6.0913599887^; - 01 
-2.1577096313^; + 00 
+1.7747816032^; + 00 
+1.6803557707^ + 00 
-7.0542734123^ - 01 



+2.1121165236^; -01 
+2.8065373289E - 01 
+2.6614220828E - 01 
+6.6514730720^ - 01 
+4.5950584109^; - 01 
+1.4752196951^; + 00 
+1.6803557707^; + 00 
+7.5925521955^ + 00 
-2.6172154724^ + 00 



-1.6655183115E - 01 
-2.1121165236E - 01 
-1.9737508953E - 01 
-4.7735386064E - 01 
-2.9366643189^; - 01 
-9.3685970173^; - 01 
-7.0542734123^; - 01 
-2.6172154724^ + 00 
+1.3630206938^ + 00 



Figures 13. String of high speed vehicles. Matrix X is 9 by 9. 
Numerical solution of stationary Riccati equation (4.10). 
The parameters = 0.1 and At = 0.1 have been used in homographic scheme. 



4.4) Control of the wave equation. 

• The fourth example is the control of the wave equation in one space dimen- 
sion 

(4.11) - = f^li{x)u,{t), xe]0,L[ 

i=i 

with homogeneous Dirichlet boundary conditions 

(4.12) y{t,0) = y{t,L) = 0. 

For solving problem (4. 11)- (4. 12), we use a spectral decomposition on the eigen- 
modes ^j{x) that are solution of the stationay problem : 

9^ $ (x) 

(4.13) -^-Q^ = , xe]0,L[ 

(4.14) $,(0) = $,(L) = 
and are classically explicited by 

(4.15) ^j{x) = V2 sin (^) , xe]0,L[ 

(4.16) A,- = ^ , j = 1, 2, 3--- 

• In practice we restrict to an approximation with the N first modes : 

N 

(4.17) y{t, x) = J2 vAt) ^j(^) . ^ > . ^ e] , L[ , 
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and with such a spectral approximation, problem (4. 11)- (4. 12) is projected with 
scalar product < • , • > 

(4.18) < u,v > = Y [ u{t)v{t)dt 

L Jo 

and the discrete formulation stands as 

(4-19) "d^ + = E ^^(^) < ^^(•) ' ^jC-) > , 1 < i < . 

i=l 

We reduce this discrete differential system to a first order one by setting 

(4.20) Yit) = en 

(4.21) C/(t) = (ni(t) , n2(t) , • • • , Umit) )^ € 

(4.22) A = c2diag(Ai, ••• , Xn) 

(4.23) Cji = < 7,(.) , >, l<^<n, l<j<7V 
and system (4.19) can be re- written as 

(4.24) lyw = (_'\ l)nt) + (0)f/(t). 

The matrices R and Q associated to the definition of the cost function J(«) 
(relation (1.6) are of order n and 2N respectively. We have chosen the follo- 
wing simple form parameterized by a = 1 and ^5 = 10 in this particular test 
case : 

(4.25) R = diag (a , • • • , a) 

(4.26) Q = diag(/3, ••• ,/5). 

• We have tested the scheme for matrices of order n — 2N — 10 and for time 
step At = 0.01 and a parameter fi = 0.001. We represent on Figures 14 to 23 
the ten different eigenvalues of Riccati matrix X{t) and for these parameters 
the scheme is stable as for as the explicit two stages Runge-Kutta scheme is 
unstable. 
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'9- 9D lis QD 19B" :3JB Sip E- -iO:- IB '120 UD 2St 



Figures 14 and 15. Wave control test. 
Two first eigenvalues of numerical solution. 




Figures 16 and 17. Wave control test. 
Third and fourth eigenvalues of numerical solution. 



fiftli eigenvalue sixth eigenvalue 




.0'. m -:j3e iff an-- o- -a- n '.Isc iff' v^-- -ais:; sek 



Figures 18 and 19. Wave control test. Fifth and sixth eigenvalues. 
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s: M:- Q \^ aEv sia 3si E '4a m -120 UD :aE 3ia aa ZHt 



Figures 20 and 21. Wave control test. Seventh and eighth eigenvalues. 




E' 40:- ai •120 ]£□ aE" 2tB .SB S!^ E - 4B Bl ' 12E ]£□ ISE" 349' :^E S(t 



Figures 22 and 23. Wave control test. Ninth and tenth eigenvalues. 



5) Conclusion. 

• We have proposed an harmonic scheme for the resolution of the matrix 
Riccati equation. The scheme is implicit, unconditionnaly stable, needs to use 
one scalar parameter and to solve a linear system of equations for each time 
step. This scheme is convergent in the scalar case. In the matrix case, harmonic 
scheme has good monotonicity properties and discrete solution tends to the 
positive solution of algebraic Ricatti equation as discrete time tends to infinity. 
We have computed first test cases of matrix square root, harmonic ocsillator, 
string of vehicles and discretized wave equation where classical explicit schemes 
fail to give a definite positive discrete solution. Our first numerical experiments 
show stability and robustness when various parameters have large variations. 
We plan to develop this work in two directions : first prove the convergence of 
the harmonic scheme in the case of Ricatti matrix equation and second construct 
a multistep version in order to achieve second order accuracy. 
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